Effect of Increasing Disorder on the Critical Behavior of a Coulomb System 
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We have performed a Monte Carlo study of a classical three dimensional Coulomb system in which 
we systematically increase the positional disorder. We start from a completely ordered system and 
gradually transition to a Coulomb glass. The phase transition as a function of temperature is second 
order for all values of disorder. We use finite size scaling to determine the transition temperature 
Tc and the critical exponent v. We find that Tc decreases and that v increases with increasing 
disorder. We also observe changes in the specific heat, the single particle density of states, and the 
staggered occupation as a function of disorder and temperature. 
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I. INTRODUCTION 

Electrons with long range Coulomb interactions dis- 
play a rich and complex behavior. In doped semi- 
conductors and disordered metals, electrons are in the 
presence of quenched disorder, and the competition be- 
tween Coulomb interactions and disorder produces a 
Coulomb glass which is an amorphous insulator. A 
great deal of effort has been expended in studying vari- 
ous thermodynamic properties of Coulomb glasses such 
as the specific heali^, the presence of a Coulomb glass 
phase transition in which the electrons are frozen into a 
highly disordered arrangemen1j2i!i£i£i£, and the Coulomb 
ga p 8 ' 9 ' 10 . Coulomb interactions between localized elec- 
trons result in the so-called Coulomb gap in the sin- 
gle particle density of states that is centered at the 
Fermi energy. Simulations have found a Coulomb gap 
in the density of states&iiiiSiiiiiiii^, and experimental 
evidence for a Coulomb gap has been seen in tunneling 
measurementsiSiiLi2iiS. 

Many of the theoretical studies of Coulomb glasses 
have been as a function of temperature. In this paper 
we will study what happens as we vary the amount of 
disorder as well as the temperature. We will start with 
an ordered system and study the effect of gradually intro- 
ducing disorder into a three dimensional system of elec- 
trons with long range Coulomb interactions. The system 
is discrete in the sense that the electrons sit on half of the 
available sites. In the ordered case the sites form a cu- 
bic lattice. The disorder is introduced in the positions of 
the sites and their deviation from the positions in a cubic 
lattice. For all values of disorder, the system undergoes a 
second order phase transition as the temperature is low- 
ered. We will study the effects on the thermodynamics 
of this phase transition as a function of disorder. 

Discretizing our Coulomb system means that it corre- 
sponds to an Ising system with long range interactions. 
An site occupied with an electron corresponds to spin- 
up and an empty site corresponds to spin-down. This is 
a very general model, and as a result relevant work has 
been done in other fields motivated by somewhat different 
physical systems. In particular there is the Ising model 
with long range interactions. Also the ordered case is 



related to work that has been done on ionic fluids near 
criticality. It is worth briefly reviewing the work that has 
been done in those fields. 

In the case of translational invariance, ionic fluids near 
criticality have been a subject of both experimental and 
theoretical investigations22i2iiS2iS. As in the case of elec- 
trons, this system is somewhat simplified by discretizing 
the system and only allowing the charges to sit on spec- 
ified sites. For ionic fluids this is known as the lattice 
restricted primitive model (LRPM)2ii22iS where there 
are equal numbers of positive and negative ions with the 
same diameter sitting on lattice sites. In the LRPM there 
is no quenched disorder. There are positive sites, neg- 
ative sites, and neutral sites (empty sites) correspond- 
ing to an Ising spin-1 model with Coulomb interactions. 
The phase diagram in the density-temperature plane has 
a second order transition line from a high temperature 
paramagnetic phase to a low temperature antiferromag- 
netic phase2ii22i2ii. This transition is in the Ising uni- 
versality class with critical exponent v — 0.63. At even 
lower temperatures there is a first order phase transi- 
tion in which the system undergoes a phase separation 
into a high density ordered phase and a low density dis- 
ordered phase. If there are no neutral sites (ionic den- 
sity p = 1), which corresponds to the antiferromagnetic 
spin-1/2 Ising model, then there is just the second order 
transition from the high temperature disordered phase to 
the low temperature ordered antiferromagnetic phase in 
three dimensions. For the purposes of this paper we are 
interested in this case where there are no neutral sites. 
Every positively charged site has a positive ion or missing 
electron, and every negatively charged site has a negative 
ion or an electron. The fact that the ionic system has a 
second order phase transition to an ordered antiferroelec- 
tric arrangement of ions2ii22i2i means that we expect the 
analogous transition to occur for the case where the elec- 
trons can sit on alternate lattice sites with no quenched 
disorder. 

Comparing the ordered and disordered extremes re- 
veals similarities and differences. Both systems undergo a 
phase transition when the temperature is lowered. In the 
ordered case, the transition is to an ordered arrangement 
of electrons occupying every other site whereas in the dis- 
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ordered case the electrons are frozen into the highly dis- 
ordered arrangement of a Coulomb glass 3 ". Both systems 
at low temperatures have a gap in their single particle 
density of states. 



As we mentioned earlier, systems with either positively 
charged sites (missing electron) or negatively charged 
sites (electron present) can be mapped onto Ising spin- 
1/2 systems. A great deal of work has been done on 
Ising models with long range interactions. The ferromag- 
netic or attractive Ising model with power law interac- 
tions that fall off as 1 /r v without quenched disorder has 
been studiedSiSSi as a function of the dimension d and 
the exponent 77 for 77 > d. However, in this paper we will 
focus on interacting electrons and so we are interested in 
the antiferromagnetic Ising model. 



The presence of quenched disorder results in an Ising 
spin glass. There has been a substantial amount of nu- 
merical effort to understand the energy of domain walls at 

T = 0,26,27,28,29,30,31,32,33 The energy Q f the domain waU 

goes as L 6 where L is the system size and the exponent 9 
is positive for systems with nonzero transition tempera- 
tures. Work on the Ising spin glass with power law inter- 
actions has been summarized in a couple of papersS^i. 
The system has a rich phase diagram in the d — 77 plane 
which can be found in Ref. 26 , where d is the dimension 
and 77 is the exponent of the power law interaction l/r n . 
The smaller 77 is, the longer the range of the interaction. 
If the range is long enough or if the dimension is large 
enough, then there is a second order phase transition 
with a transition temperature Tc > 0. The critical ex- 
ponents are different in the long range and short range 
regimes. The exponent 9 depends continuously on 77 in 
the long-range region (9 = d — 77), and is independent of 
7/ in the short-range regime 26 . This indicates that the 
critical exponents also depend continuously on 77 in the 
long-range region, and are independent of 77 in the short- 
range regime 2 ^. Katzgraber and Young have done Monte 
Carlo simulations of an Ising spin glass in one dimension 
with long range interactions 2 ^ 7 .. They chose a value for 
77 where the system has a second order spin glass transi- 
tion, and they find that v — 10/3. 



In this paper we will be concerned with what happens 
to thermodynamic quantities as we systematically intro- 
duce disorder into a three dimensional system of electrons 
with long range Coulomb interactions. The disorder is in- 
troduced into the placement of sites where the electrons 
can sit. The paper is organized as follows. In section 
II we present the Hamiltonian and describe our Monte 
Carlo simulation. In section III we present the quantities 
that we measure. In section IV we present our results, 
and we give our conclusions in section V. 



II. CALCULATION 

A. Hamiltonian 

Let us start by considering the completely disordered 
case which is known as a Coulomb glass. The essential 
physics of the Coulomb glass is the presence of both dis- 
order and long range Coulomb interactions between elec- 
trons. The Hamiltonian often studied for the Coulomb 
glass is^ii 

g = E^+E ^" x)( " J " J " ) (i) 

where we set the charge e = 1, rii = ±1 is the number 
operator for site i, (pi is the onsite energy, ry = 1 7^ — fj\, 
and K is a compensating background charge making the 
whole system charge neutral. Such a Hamiltonian de- 
scribes a lightly doped semiconductor, in which the im- 
purity sites are far enough apart that the overlap between 
sites can be neglected. In most of the early work on the 
Coulomb glass (e.g., Refsi 4 i 6 i 35 ), the sites are chosen to 
form a periodic lattice, and the disorder is present in the 
form of random onsite energies. For an ordered system, 
the onsite energy <pi is a constant. One could imagine 
gradually introducing disorder by allowing fa to be cho- 
sen from a distribution whose width gradually increases. 

However, the presence of random onsite energies makes 
numerical analysis difficult, since even in the high tem- 
perature state the average occupation of a site is not 
zero. This makes the search for a phase transition diffi- 
cult; there is no obvious order parameter which becomes 
nonzero at the transition. For our numerical analysis, it 
is more convenient to take the disorder to be entirely in 
the location of the sites. This changes the symmetry of 
the Hamiltonian from having onsite disorder to having 
disorder in the interaction between sites because the dis- 
tance between sites varies. For many quantities these two 
models give similar results. For example studies of the 
specific heat in Coulomb glasses have compared having 
disorder in the onsite energy to having a completely ran- 
dom displacement of site o 2 i 36 . They find that both mod- 
els produce qualitatively similar results with some quan- 
titative differences. However, the existence and nature of 
the phase transitions is different in the two models^L 3 ^. 
In particular there is always a phase transition no mat- 
ter how wide the distribution of the site placement is, 
whereas there is no phase transition if the width of the 
distribution of the onsite energy fa is larger than a critical 
valued. Mobius 3 ^ has argued that such a critical value 
must exist, even if it is vanishingly small, since there 
is a phase transition when there is no disorden4°« while 
there is no clear evidence for a transition when there is 
substantial onsite disorder. This implies that long range 
order is destroyed by both onsite disorder and thermal 
fluctuations. 

A number of previous simulations have used the form 
for the Hamiltonian with disorder in the placement of the 
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site a 2 i 3 i 36 i 41 . In the case of half filling there is a particle- 
hole symmetry, and the phase transition is associated 
with the development of a nonzero Edwards-Anderson 
order parameter—. We therefore rewrite the Hamiltonian 
(taking K = 1/2) to look like that of a spin glass, 

s,s i (2) 



H 



i>j 



Si — 1 (—1) will denote an occupied (unoccupied) site. 

We have simulated three dimensional systems of linear 
size L = 4, 6, and 8. We place N = L 3 sites in the sys- 
tem. We have only considered the case of half filling in 
order to take advantage of the spin-flip symmetry. For 
the ordered case the sites form a cubic lattice. In the 
ground state, every other site is occupied; the occupied 
sites form a face centered cubic (FCC) lattice. We can 
gradually introduce disorder by allowing the deviation of 
a site from its position in a cubic lattice to be chosen 
from a Gaussian distribution with a standard deviation 
of a. This gives the radial distance from the cubic lattice 
site. The angular coordinates of the site are chosen ran- 
domly using a uniform distribution. The ordered case 
corresponds to a — 0. a = 1 corresponds to the very 
disordered case with a standard deviation equal to the 
cubic lattice constant a. We also considered completely 
random arrangements of sites where the x, y, and z co- 
ordinates of each are chosen from a uniform distribution. 
We call this the "uniform random" case. We find no qual- 
itative difference and only a slight quantitative difference 
between the uniform random case and the a = 1 case in 
quantities such as the single particle density of states, the 
specific heat versus temperature, and the Binder's g. So 
we will not make much mention of the uniform random 
case. 

We use infinite periodic boundary conditions in which 
the simulation box is infinitely replicated in all directions 
to form a lattice. As a result, an electron on a given site 
interacts with other electrons and all their images via the 
Coulomb interaction. To handle this, we use an Ewald 
summation technique^ which replaces the Hamiltonian 
in Eq. (J2J) with the following effective interaction between 
sites: 

l<i<j<N v ' i=l 

where L is the linear size of the simulation box, N is the 
number of sites, the charge qi — Si/2, and the function 
i/j(r) is given by 

erfc(a|r + n|) 



m = £ 



\r + ri\ 



^£^ e H 27 



in which 



2 f x _ 2 
erfc(x) = 1 -= / e 1 dt 

V 71 " Jo 



(4) 



(5) 



is the complementary error function and 
"erfc(a|n|) , 1 „-^ n *ic? 



Note that 



2a 

7^ 



A = lim 

|r|-»0 



(6) 



(7) 



The sum over n in Eq. JBJ is a sum over all simple cu- 
bic lattice points with integer coordinates n = (l,m,n). 
These are the coordinates of the images of the simulation 
box. The parameter a is a convergence factor that is ad- 
justed to maximize the rate of convergence of the sum. 
We have omitted a positive term in the Hamiltonian © 
that is proportional to the square of the net dipole mo- 
ment of the configuration^. This omission is equivalent 
to the boundary condition in which the infinite sphere 
of our system and its images is surrounded by a perfect 
conducting medium. We have done some runs with the 
dipole term and find no qualitative difference and only a 
slight quantitative difference compared to the case with 
no dipole term. So in this paper we will present the result 
of simulation runs which omit the dipole term. 



B. Monte Carlo Simulation 

We have used a Monte Carlo heat bath algorithm. We 
keep a table of the potential energy at each site. Each 
electron is looked at sequentially and moved to one of the 
available N/2 + 1 sites (its own site or one of the available 
N/2 unoccupied sites), chosen with a Boltzmann proba- 
bility. If the site chosen is the electron's original location, 
the potential energies are unchanged; if the electron hops 
to a new site, we update all the potential energies. If the 
electron chooses its initial site, which it does with high 
probability at low temperatures, we do not have to re- 
compute the potential energies. This speeds up the simu- 
lation considerably, partially compensating for the much 
longer equilibration times needed at low temperatures. 
Our longest run (for L — 4 at T — 0.01 and a = 0.5) 
had 3 x 10 6 Monte Carlo steps per electron. Depending 
on the system size and temperature, the sample averages 
involved between 5 and 190 disorder configurations. 



III. MEASURED QUANTITIES 

A. Binder's g and equilibration criteria 

The Edwards-Anderson order parameter alluded to 
above quantifies the extent to which spins or site oc- 
cupations are frozen. It is defined as q = [(Si) 2 ] ; we 
will denote thermal averages by ( ... ) and disorder av- 
erages by [ ... ]. Thermal averages sum over fluctuations 
in the positions of the electrons weighted with the cor- 
rect Boltzmann probability; disorder averages sum over 
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different arrangements of the sites. We can see from the 
definition of q that if the spins are frozen, then the av- 
erage orientation of a spin will have a nonzero thermal 
average and q will be finite. This is why q can be thought 
of as the order parameter of the phase transition. 

We can generalize the Edwards-Anderson order pa- 
rameter to a finite time overlap in either of two ways 4 ^. 
The first way computes the overlap between two replicas: 



(8) 



where the superscripts refer to different replicas. The 
two replicas are identical in their disorder, i.e., the place- 
ment of the sites, but differ in the initial positions of the 
electrons. The other way uses the same replica at two 
different times 



(9) 



If the time difference r is sufficiently large that the elec- 
tron configurations at t and t + r are essentially uncor- 
rected, q t (t,T) will give the same result as the replica 
overlap. 

We use the moments of the overlap to define Binder's 
g which is a parameter that is related to the phase tran- 
sition. First we define g(r) by^Mi 



where 



2 I [(9 2 W)] 2 , 



(? n (r)> = -£«"(*) 



(10) 



(11) 



We will use gt (g r ) to denote the result of using q t (q r ) 
Binder's g is given by 



= lim gir) 



which we will approximate by 



9(r) 



(12) 



(13) 



for some measurement time r large enough that the con- 
figurations are essentially uncorrelated so that gt and g r 
agree. We have used this fact to monitor equilibration by 
simulating two replicas that have the same placement of 
sites but different spin configurations. 43 Typically g r (t) 
increases with time to the equilibrium value, whereas 
gt(t) decreases to the equilibrium value. The two meth- 
ods agree when the system has reached equilibrium. Our 
criterion for equilibration was that the values of g r and gt 
agreed to within 0.1. We only present results for systems 
which meet this criterion. We should caution that this 
criterion can be met even though the system may still be 
slowly aging. 



Binder's g provides a way to monitor the phase tran- 
sition. At high temperatures, the distribution of q tends 
to a Gaussian so that g — > 0, whereas the order pa- 
rameter, and hence <?, become nonzero as the temper- 
ature approaches the phase transition temperature Ta- 
li we make the assumption of one parameter scaling, 
then the only relevant length is the correlation length 
£ ~ (T — Tc)~ u where v is the critical exponent associ- 
ated with £. So all lengths, including L, can be scaled 
by £. Since g is dimensionless, we expect that it should 
satisfy a scaling form 43,44 



g(L,T)=g(L l /»(T-T c )). (14) 

Thus at the critical temperature, g(L,Tc) should have 
the same value independent of the system size L (as 
long as L is sufficiently large for finite size scaling to 
apply) i 3 - 4 ! 



B. Specific Heat 

There are two ways to calculate the specific heat 
Cv{T). The first way uses the variance of the energy 
fluctuations: 



1 



Nk R T 2 



[(E 2 ) {Ef 



(15) 



where E is the average energy per electron, N is the 
number of electrons, and ks is Boltzmann's constant. 
The other way to calculate the specific heat is to take the 
derivative of the energy with respect to temperature. We 
can approximate the derivative by a finite temperature 
difference 



C v (Ti) 



d[(E)} 



dT 



i(E(T +1 ))} - [(E(Ti))] 



T, 



(16) 



We found that the specific heat calculated in these two 
ways agreed quite well. Notice however, that if the slope 
of E versus T is increasing as temperature increases, 
then the specific heat calculated by the finite difference 
method will underestimate Cy which is actually the slope 
of the tangent to the energy curve. Similarly, if the slope 
of E versus T is decreasing as the temperature increases, 
the finite difference method will overestimate Cy. 



C. Staggered Occupation 

Since an unoccupied site on a cubic lattice corresponds 
to a down Ising spin and an occupied site to an up Ising 
spin, the FCC crystalline phase corresponds to a maxi- 
mum in the magnitude of the staggered occupation M s . 
The staggered occupation is defined by£ 



M s = i^(-l) i+ ^5.. 



i+j+k 



(17) 
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where i, j, k are the integer coordinates of the sites in a 
cubic lattice in units of the lattice constant a. So a site 
coordinate (a;, y, z) = (ia,ja, ka). Since disorder is intro- 
duced through a distribution in the position of the sites 
with respect to the cubic lattice sites, we can still use eq. 
(|17J) to calculate the staggered occupation in the pres- 
ence of disorder by regarding i, j, and k as coordinates 
of the center of the unit cell where the site is located. 
It is useful to plot the staggered occupation distribution 
P(M S ) versus M s in order to see the extent of the "crys- 
talline" order. In order to compare different system sizes, 
we normalize the staggered occupation to range from —1 
to +1, and the area under the curve is normalized to 1. 



D. Single Particle Density of States 

We have calculated the single particle density of states 
N(E) at various temperatures. N(E) is the distribution 
of potential energies at single sites due to interactions 
with all the other sites. In other words, we can write the 
Hamiltonian in the form of an Ising model: 



1.0 & 



H = 



2 E JijSiSj 



(18) 



where £7, is the single site energy or "local field" and is 
given by 



Ei — - ^2 Jij Sj 



(19) 



N(E) is the thermally averaged and disorder averaged 
distribution of Ei. 



IV. RESULTS 

A. Second order melting transition in the ordered 

case (a — 0) 

We consider the case where the translational degrees 
of freedom are discrete and the electrons can only sit on 
designated sites. This is equivalent to a long range Ising 
model in three dimensions. If the sites are ordered and 
are lattice sites, it is also equivalent to the lattice re- 
stricted primitive model (LRPM) in the completely filled 
case (ionic density p = 1)2^. We find that discretiza- 
tion produces a second order phase transition regardless 
of the amount of positional disorder. In this section we 
present evidence that the ordered case with a = under- 
goes a second order crystallization transition to an FCC 
lattice as the temperature is lowered. This result is con- 
sistent with the second order transition found by Mobius 
and Rofiler— from numerical simulations of a half-filled 
system on a cubic lattice with Coulomb interactions. It 




Temperature 



FIG. 1: (a) g r (L,T) vs. T for a = 0. (gt(L,T) vs. T is 
virtually identical.) The data for L — 4 is averaged over 
190 runs, L = 6 is averaged over 67 runs, and L = 8 is 
averaged over 45 runs. The solid lines are guides to the eye. 
(b) g(L, T) for a = scaled using g (l 1/v (T-T c j) with 
T c = 0.128 ± 0.001 and v = 0.55 ± 0.1. 



is also consistent with the second order transition in the 
LRPM model for the fully occupied case (p = 1)2344. 

In Figure ^ we show g versus the temperature T at 
L = 4, 6, and 8. The point where these curves cross 
yields a transition temperature of Tc = 0.127. Notice 
that above the transition g dips down and acquires neg- 
ative values. This behavior has been seen in the case of 
a 3-state ferromagnetic Potts model in three dimensions 
which undergoes a first order phase transition 4 - . How- 
ever, in that case the value of g at the minimum scaled 
as g(T m i n ) ~ — L d , whereas in our case g(T m i n ) appears 
to saturate at large L. The negative values of g can re- 
sult if the distribution P(q) is nongaussian with finite 
weight at q ^ corresponding to long lived occupations 
of some sites. A very simple delta function distribution 
that illustrates this is 



P(q)=a 5(q) 



1 



1 - a a 



5(q-a )+ — ^—2- S(q + a ) 

(20) 

where a is a parameter with values between and 1, 
and a is a constant. For 2/3 < a < 1, this distribution 
yields g < 0. 

We initially thought that the transition might be first 
order. One of the signatures of a first order melting tran- 
sition is coexistence of the liquid and crystalline phases 
at the melting temperature. We looked for evidence of 
coexistence by examining the distribution P(M S ) of the 
staggered occupation. Coexistence would produce three 
peaks in P(M S ) versus M s : a central peak and two side 
peaks symmetrically placed with respect to M s = 0. The 
central peak corresponds to the high temperature liquid 
phase and the side peaks correspond to the FCC crys- 
talline phase. Furthermore, at the transition tempera- 
ture for a first order transition, the three peaks would 
become narrower and higher with increasing system size. 
On the other hand, if the system is cooled through a sec- 
ond order transition, the high temperature central peak 
in P(M S ) is replaced by two peaks symmetrically placed 
about M s = 0. These peaks do not become sharper 
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FIG. 2: The distribution P(Ms) of the staggered occupation 
Ms for L = 8 at a = at various temperatures. The central 
peak is highest at T = 1 and gradually decreases as T de- 
creases. The two side peaks begin to appear in the vicinity of 
Tc and become more pronounced as T drops below Tc ■ There 
is no temperature where 3 peaks are present, indicating that 
the transition is not first order. The data was the result of 
averaging over 35 runs. 




Temperature 

FIG. 3: g r (L = 8,T) vs. T for a = (45 runs), a = 0.1 (10 
runs), a = 0.2 (5 runs), a — 0.3 (15 runs), a — 0.4 (115 runs), 
a = 0.5 (45 runs), and a = 1 (108 runs). (g t (L = 8,T) vs. 
T is virtually identical.) The number of runs in parentheses 
is the number of runs that were averaged to obtain the data. 
The solid lines are guides to the eye. 



ture. 



with increasing system size, but the width of the dis- 
tribution is expected to decrease with increasing sys- 
tem size as L~$l v where f) and v are the critical expo- 
nents defined by M s ~ \T - T c f and £ ~ (T - T c y" ' 47 
Figure [3 shows the distribution P(M S ) of the staggered 
occupation at various temperatures. Notice that in the 
vicinity of the melting temperature there are only two 
symmetrical side peaks. This implies that the transi- 
tion is a second order phase transition. Furthermore we 
find that the value M s ^ max , where P(M S ) has a maxi- 
mum, decreases with increasing system size at Tc and 
goes as M s ^ max ~ L~ 0,6 . This is also consistent with 
a second order transition. In the vicinity of the phase 
transition where P(M S ) has 2 peaks, we can define the 
width M StW idth of the distribution as the nonzero value 
of M s where P{M s<width ) = P(M S = 0). We find at T c 
that M SjW idth is linear in L and can be fit to the form 
M SiW idth — A — mL where A and m are constants that 
are temperature dependent. At T — 0.128 which is close 
to T c , A = 1.1 and m = 0.027. Notice that M StWidth 
does not appear to follow the form M StWi( ith ~ L~P/ V ', 
but we would need more than 3 values of L to accurately 
determine if there is a discrepancy with the scaling form. 

First order phase transitions are often characterized by 
hysteresis upon heating and cooling. We have looked for 
hysteresis by cooling and then heating the system, and 
examining the resulting curves of g versus T as well as the 
specific heat Cy versus T. We find no hysteresis which 
is further evidence against a first order phase transition. 

To summarize, the ordered case (a = 0) undergoes a 
second order phase transition as a function of tempera- 



B. Critical Behavior 

We have determined the critical exponent v and the 
transition temperature Tc as a function of the disorder 
a through the finite size scaling of g(L, T)m^ In Figure 
E]we plot g(L = 8,T) versus T for various values of a. 
Notice that the transition region moves to lower tempera- 
tures with increasing disorder. This reflects the decrease 
in Tc with increasing a. The transition temperature cor- 
responds to the temperature where the curves of g{L, T) 
versus T for all sizes cross. Examples are shown in Fig- 
ure^and Figure^ To more accurately determine Tc, we 
use the scaling hypothesis to collapse the data for a given 
value of a onto a single curve as shown in Figure Tc 
and v are used as adjustable parameters to collapse the 
data. The values of v and Tc at various values of a are 
given in table||| We can estimate the errors in the critical 
temperature and the critical exponent v by how well the 
curves can be made to collapse. The errors given in the 
table also include our estimate of the effects of aging. In 
other words, the error bars include our estimate of how 
the values might change if we were to run longer at low 
temperatures or cool more slowly. In Figures |3] and we 
plot Tc and v versus a. 

We can see that v increases from v = 0.55±0.1 at a = 
to v = 1.30 ±0.2 at a = 1. The value of v in the ordered 
case (a = 0) lies between the classical value (v = 0.5) 
and the value for the ordered short ranged Ising model 
{v = 0.63) 48 . Within the error bars, our value is consis- 
tent with both universality classes, and therefore cannot 
differentiate between them. Mobius and Rofiler 40 studied 
a half-filled system on a cubic lattice with Coulomb in- 
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FIG. 4: (a)-(c) g r {L,T) versus T for a = 0.3, 0.4, and 1.0 
at L = 4, 6, and 8. The solid lines are guides to the eye. 
(g r (L — 8, T) vs. T is virtually identical.) The number of 
runs in parentheses is the number of runs that were aver- 
aged to obtain the data, (d) g(L,T) for a = 0.3 scaled using 

g (T - T c )) with T c = 0.085±0.005 and v = 0.71±0.1. 

(e) g(L,T) for a = 0.4 scaled using g ^L 1/v (T - T c )\ with 
T c = 0.045 ± 0.01 and v = 1.05 ± 0.1. (f) g(L, T) for cr = 1 
scaled using g [l 1/v (T - Tcfj with T c = 0.028 ± 0.01 and 
v= 1.30 ±0.2. 



TABLE I: The values of Tc and v for different valuse of a. 



a 


T C 


V 


0.0 


0.128 ±0.005 


0.55 ±0.1 


0.1 


0.123 ±0.005 


0.57 ±0.1 


0.2 


0.110 ±0.005 


0.61 ±0.1 


0.3 


0.085 ± 0.005 


0.71 ±0.1 


0.4 


0.045 ± 0.01 


1.05 ±0.2 


0.5 


0.030 ± 0.01 


1.35 ± 0.2 


1.0 


0.028 ± 0.01 


1.30 ± 0.2 



teractions and found v = 0.635(10) which agrees with 
the value for the Ising model. Our simulations differ 
from those of Mobius and Rofiler in that we used the 
Ewald summation to take into account the fact that the 
Coulomb interaction extends beyond the size of the sys- 
tem while they did not. As we mentioned earlier, the 
order-disorder transition in the LRPM model on a sim- 
ple cubic lattice belongs to the Ising universality class^. 




a 

FIG. 5: The critical exponent v versus the disorder a. The 
solid line is a guide to the eye. 



The completely filled LRPM model with ionic density 
p = 1 is equivalent to our a = case. In addition Lui- 
jten et al. did Monte Carlo studies of the restricted primi- 
tive model (RPM) which has equal numbers of oppositely 
charged ions with equal diameters and with Coulomb in- 
teractions in three dimensions 4 ^. These grand canonical 
simulations of the RPM used a finely discretized lattice 
where the ionic diameters were 5 times larger than the 
lattice spacing, and they found the Ising of the critical 
exponent v = 0.63(3). 

In the disordered case (a = 1) our value for v = 
1.3 ± 0.2 differs from the value of v = 0.75+o^ obtained 
earlier' 5 . Again this is probably due to the fact that we 
used Ewald summation whereas the previous work did 
not. In addition we were able to do longer runs at low 
temperatures than the previous work. 

The transition temperature decreases from Tc = 
0.128 ± 0.005 at a = to T c = 0.028 ± 0.01 at a = 1. 
The value of Tc — 0.128 at a = is consistent with the 
temperature of the peak in the specific heat found pre- 
viously by Mobius and Rofiler 40 . Within the error the 
value of Tc = 0.028 ± 0.01 at a = 1 is consistent with 
the previous value of T c = 0.043±g;gg| found by Grannan 
and YvA 

It is interesting that Tc is much lower than the char- 
acteristic energies of the system which are of order unity. 
This is especially true for large values of the disorder. 
The reason for this was given by Grannan and Yu^ and 
is as follows. At the temperatures of our simulations, 
nearby pairs of sites will with high probability consist of 
an occupied and an unoccupied site. Since these strongly 
coupled pairs of sites are close together, they are guaran- 
teed to have small dipole moments. Therefore, they will 
interact weakly with the rest of the system, remaining 
active down to temperatures much lower than the bare 
interaction energy. 
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FIG. 6: Transition temperature Tc vs. a (o), temperature 
T m ax of the maximum of dN(E = Q)/dT vs. a (□), and the 
temperature T pea k of the maximum in the specific heat vs. a 
(0). The solid lines are guides to the eye. 
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Temperature 

FIG. 7: The specific heat Cv versus T in units of ks for 
L = 8 for a = (45 runs), a = 0.1 (10 runs), a = 0.2 (5 
runs), a = 0.3 (15 runs), a = 0.4 (95 runs), a = 0.5 (45 
runs), and a = 1 (108 runs). The number of runs averaged 
over is indicated in parentheses. The solid lines are guides to 
the eye. 



C. Specific Heat 

In Figured we plot the specific heat versus tempera- 
ture for various values of a for L = 8. We see that in the 
ordered case (a = 0) Cy exhibits a sharp peak centered 
at Tq- As the disorder increases, the peak broadens and 
eventually becomes a broad bump with a maximum at a 
temperature above Tc- For example, for a = 1, Cy has 
a maximum at T = 0.07 whereas Tc — 0.028. In Figure 
Elwe compare the temperature T pea k of the maximum in 
the specific heat with Tc for various values of a. We see 
that Tp ea k matches well with Tc for a < 0.3. For larger 
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FIG. 8: The specific heat Cy versus T in units of ks for 
L = 4, 6, and 8 for (a) a = and (b) a = 1. The specific 
heat is calculated from the variance in the energy fluctuations. 
The specific heat calculated from the derivative of the energy 
with respect to temperature is similar. The number of runs 
averaged over is indicated in parentheses. The solid lines are 
guides to the eye. 



values of a, T pea k > Tc- Spin glasses also have a max- 
imum in their specific heat at a temperature above the 
spin glass transition temperature^. For the three dimen- 
sional Coulomb glass where the disorder is in the onsite 
energy rather than in the positions of the sites, Mobius et 
al. found that as the width in the distribution of onsite 
energies increased, the temperature T pea k of the maxi- 
mum in the Cy also increased^. However, in the cases 
of onsite disorder that they considered, the maximum 
does not signify a transition since the existence of a phase 
transition in the presence of onsite disorder has not been 
established. The maximum in Cy must be present since 
Cy(T) goes to zero at the extremes T — > and T — > oo, 
implying that there must be a maximum in between these 
extreme* 39 . Furthermore, even without Coulomb inter- 
actions but with a large amount of onsite disorder, there 
would be a maximum in the specific heat consisting of a 
superposition of the Schottky specific heats of two level 
systems with randomly distributed excitation energies 3 ^. 

To show the size dependence of the specific heat, in 
Figure |H1 we plot Cy versus T for different system sizes 
at fj = and at a = 1. In the ordered case the specific 
heat peak becomes sharper as L increases while in the 
disordered case, the broad bump is only weakly depen- 
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FIG. 10: N(E) versus E for L — 8 for various values of a 
at temperatures in the vicinity of Tc- Shown are a — 
(T = 0.128, 45 runs), a = 0.1 (T = 0.122, 10 runs), cr = 0.2 
(T = 0.105, 5 runs), a = 0.3 (T = 0.085, 15 runs), a = 0.4 
(T = 0.045, 115 runs), a = 0.5 (T = 0.030, 45 runs), and 
cr = 1 (T = 0.0275, 108 runs). The temperatures and the 
number of runs averaged over is indicated in parentheses. The 
solid lines are guides to the eye. 



FIG. 9: N(E) versus E for L — 8 at various temperatures. 
The solid lines are guides to the eye. (a) a = 0. The data was 
averaged over 45 runs, (b) a = 1. The 3 lowest temperatures 
were averaged over 108 runs and the 3 highest temperatures 
were averaged over 16 runs. 

dent on system size. 

D. Single Particle Density of States N(E) 

In a Coulomb glass the long range Coulomb interac- 
tions between localized electrons produce a Coulomb gap 
in the single particle density of states that is centered 
at the Fermi energ} ' 8 ^ 10 . The Coulomb gap makes the 
ground state stable with respect to single electron hops. 
The ordered case also has a gap but for a somewhat dif- 
ferent reason. In the ground state of the ordered case 
where there is an FCC lattice, the potential energy or 
local field is the same for each occupied site. The local 
field is equal and opposite for the unoccupied sites. This 
leads to an N(E) with two delta functions symmetrically 
placed about E = 0. In finite size systems at finite tem- 
peratures these delta functions broaden into finite height 
peaks due to thermal fluctuations and the formation of 
ordered domains. 

In Figure we show the density of states N(E) for 
single particle excitations at various temperatures for a — 
and for a = 1. Because of strong electron-electron 
correlations, the density of states at zero energy starts to 



decrease at about TTc in the ordered case (a = 0) but 
at a temperature about an order of magnitude above Tc 
in the strongly disordered case (cr = 1). In Figure ITOI we 
show N(E) at or near Tc for various values of cr for L=8. 
We see that at Tc the gap appears nearly fully formed 
for cr = 1 but not for a = 0. In the ordered case the finite 
density of states at E = is possibly due to domains. 

As we can see from the figures, at finite temperatures 
the gap in the density of states is partially filled, and 
the density of states does not vanish at the Fermi energy. 
This has been seen in previous simulationsiii*i2iiiiii^. 
Tunneling measurements of the Coulomb gap have also 
seen that it fills in with increasing temperatur e) 1 ^ 1 ? . The 
exact form of N(E, T) is not known, but for strong 
disorder some have arguedii*i2iii that its low tem- 
perature asymptotic behavior is described by N(E = 
0,T) ~ T d ~ x . However, some simulations^ have found 
a stronger temperature dependence, i.e., N(E = 0, T) ~ 
T A with A > (d— 1). For d= 2, Sarvestani et al~ found 
A = 1.75 ± 0.1, and for d = 3, A = 2.7 ± 0.1. 

In Figure ITT1 we show our results in a log-log plot of 
N(E — 0, T) for various values of a at L = 8. 

At low temperatures the curves are quite straight on a 
log-log plot. So we can fit the low temperature part of 
these curves to a power law form N(E — 0, T) ~ T A . The 
fits are shown as solid lines in Figure 1111 and in Figure 
[TI\ We plot A as a function of a in Figure ^| We find 
that A varies between 3 to 16 and is always greater than 
cZ — 1 = 2 since d = 3. Even in the case of uniform 
disorder (uniform random), A = 4.8. The large value 
of A for cr = is not entirely surprising since Figure 



10 



o 

ii 

LU 



10' 
10° 
10" 2 
1CT 4 

10- 6 

10" r 




o & 3 



O o=0.0 
□ o=0.1 
O o=0.2 
A o=0.3 
< o=0.4 
V o=0.5 
> o=1.0 



10" 



10" 1 

Temperature 



10" 



FIG. 11: Log-log plot of N(E = 0) versus T for L = 8 for 
(7 = (45 runs), a = 0.1 (15 runs), a = 0.2 (10 runs), a = 0.3 
(10 runs), a = 0.4 (115 runs), a = 0.5 (45 runs), and a — 1 
(108 runs). The number of runs averaged over is indicated in 
parentheses. The solid lines are power law fits to the data at 
low temperatures. 
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FIG. 12: Log-log plot of N(E = 0) versus T at low tempera- 
ture for L = 8 for a — (45 runs), a = 0.1 (15 runs), a = 0.2 
(10 runs), a = 0.3 (10 runs), a = 0.4 (115 runs), a = 0.5 (45 
runs), and a = 1 (108 runs). The number of runs averaged 
over is indicated in parentheses. The solid lines are power law 
fits to the form N(E = 0, T) ~ T A . 



fTUl shows that N(E = 0, T) is larger for a = than 
for any other value of the disorder in the vicinity of the 
transition temperature. Since N(E — 0, T) goes to zero 
as the temperature goes to zero, N(E — 0, T) has the 
farthest to go for a = 0. Even though Tq is largest for 
(7 = 0, the ratio N(E = 0,T)/Tc is largest for the case 
of no disorder, and so it is consistent that the exponent 
A is largest for the case of no disorder. 

We plot the data from Fig. ^] on a linear plot in Fig- 
ure El where we see S-shaped curves. We can see that 
N(E — 0, T) rises much more steeply for small values of 




FIG. 13: The power A versus a for L 
guide to the eye. 



The solid line is a 




0.2 

Temperature 
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FIG. 14: Linear plot of N(E = 0) versus T for L = 8. Data 
is the same as is shown in Fig. 1111 The solid lines are guides 
to the eye. 



disorder than for large values of disorder. The steepest 
part rise for the ordered cases (a < 0.3) occurs approxi- 
mately at Tc- We can quantify this by taking a deriva- 
tive dN(E = 0)/dT that can be approximated by a finite 
difference: 



dN(E = 0) 



dT 



N l+1 (E = 0) - Ni(E = 0) 



T=Ti 



Ti 



i+i 



T, 



(21) 



The result is shown in Figure where we compare Tc 
with the temperature T max where dN(E — 0)/dT is a 
maximum. We see that T max follows Tc for a < 0.3 but 
lies above Tc for larger values of the disorder a. 

Efros and Shklovskii2iia argued that at T = the 
Coulomb gap in the single particle density of states of a 
fully disordered system should scale as N(E) ~ \E — Ep\ s 
where 8 = d — 1 and d is the dimension of the system. 
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FIG. 15: Log-log plot of N(E) versus \E - E F \ for L = 8 
for a = (T = 0.120, 45 runs), cr = 0.1 (T = 0.111, 15 
runs), a = 0.2 (T = 0.08, 20 runs), a = 0.3 (T = 0.055, 15 
runs), a = 0.4 (T = 0.0335, 115 runs), a = 0.5 (T = 0.0275, 
45 runs), cr = 1.0 (T = 0.0300, 108 runs). The temperatues 
are all below Tc- The solid lines are fits to a power law 
[N(E) - N(0)] ~\E- E F \ S for values of E very close to the 
Fermi energy Ef, i.e., \E — Ef\ < 0.04. The plots include 
N(E) values for E above and below E F . 



Some subsequent workiiiisiS has supported this form for 
the density of states, though some simulations^iiSii^i 
have found a steeper energy dependence, i.e., d > 

5 > (d — 1) in two 5 i 12 i 15 i 51 and threo 12 ! 15 ' 51 dimensions. 
Efros 52 included two-electron transitions in calculating 
the density of states of a Coulomb glass and proposed 
the exponential form N(E) ~ exp [-\E /(E - Ep)] 1 / 2 ] 
where E Q is a constant. The physical reason for such a 
sharp gap is the formation of polarons in which an oc- 
cupied site tends to have unoccupied sites nearby and 
vice-versa. Some simulations^ have found support for 
this exponential form, while otheraS^i have not. 

According to the theorji^iifi, N(E) ~ \E - Ep^" 1 
in the limit E — > Ep. In Figure IT51 we plot our data 
for [N(E) - N(E F )} versus \E - E F \ on a log-log plot. 
(Since we are at finite temperatures, we have subtracted 
off N(Ep).) We fit the low energy data in the vicin- 
ity of the Fermi energy Ep to a power law of the form 
N(E) ~ \E—Ep\ 5 for various values of a at temperatures 
below Tc. For the case of strong disorder (cr = 1), we 
find 6 = 2.65 ± 0.2 which agrees with the previous values 
of 5 = 2.6±0.2 found by Mobius et al. 51 and 5 = 2.7±0.1 
found by Sarvestani et al. 15 . It disagrees with the value 
oi 5 = d — 1 = 2 predicted by Efros and Shklovski&iS 
and with the value S = 2.38 found by Li and Phillips 12 . 
In the case of no disorder, the curvature is very close to 
quadratic and we find 6 = 2.1 ± 0.2. In Figure IT^l we 
plot the exponent S versus the disorder a. We see that 

6 increases and then saturates with increasing disorder. 
The estimated error of ±0.2 in 6 does not come from the 
fit to the data, so much as from the fact that the finite 
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FIG. 16: The exponent 8 versus a from the fits to the power 
law [N(E) - N(0)] ~\E- E F \ 6 in Fig.EHl The solid line is a 
guide to the eye. The error in S is approximately ±0.2 due to 
finite size and finite temperature effects described in the text. 
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FIG. 17: Log-log plot of -l/ln[N(E) - N(E F )] versus \E - 
E F \ for L = 8 for a = (T = 0.120, 45 runs), a = 0.1 
(T = 0.111, 15 runs), a = 0.2 (T = 0.08, 20 runs), a = 0.3 
(T = 0.055, 15 runs), a = 0.4 (T = 0.0335, 115 runs), a = 0.5 
(T = 0.0275, 45 runs), a = 1.0 (T = 0.0300, 108 runs). 
The temperatues are all below Tc- The solid lines are fits to 
-l/ln[N(E)-N(E F )] ~ \E-E F \"' for values of E very close 
to the Fermi energy Ef, i.e., \E — E F \ < 0.04. The slope of 
each line gives 7. The plots include N(E) values for E above 
and below Ef. 



temperature affects low energies E ~ kT. There are also 
finite size effects^ that affect low energies E ~ 1/2L, 
though finite size effects for L > 6 are quite small (less 
than 1%). 

We have checked to see if our data provides evidence for 
the exponential form N{E) - exp [-\E /(E - Ep)] 1 / 2 ] 
proposed by Efros^. In Figure El we show a log- 
log plot of -l/ln[JV(J5) - N(E F )} versus \E - E F \ for 
various values of a. (Since we are at finite temper- 



12 



0.70 

?- 

I 0.60 

Q. 

C 

CD 

§ 0.50 

Q. 

X 
CD 



0.30 




a 

FIG. 18: The exponent 7 versus a from the fits to the ex- 
ponential form [N(E) - N(E F )] ~ exp \-E a /(E - £f)| 7 in 
Fig- El The solid line is a guide to the eye. 

atures, we have subtracted off N(Ep).) If N(E) ~ 
exp \—\E / (E — Ep)] 1 ^ 2 ] were a good description of the 
density of states, then the curves in Figure El would be 
straight lines with slopes of 1/2. Since the exponential 
form presumably only describes the density of states in 
the vicinity of Ep, we have fit lines through the points 
corresponding to \E\ < 0.04 assuming the more general 
form [N(E) - N(E F )} ~ exp - \E a /(E - E F )\~' . The 
slope of the lines in Figure El correspond to the expo- 
nent 7. We plot 7 versus a in Figure El The values of 
7 fluctuate around 1/2, but the large curvature of the 
trajectories in Figure 1171 do not lend strong support to 
the exponential form of the density of states. 

Analytical theoriesiii4 of the Coulomb glass predict 
that the finite temperature density of states N(E = 0, T) 
at the Fermi energy (Ep = 0) should be proportional to 
the zero temperature density of states N(E, T — 0) at an 
energy E = k B T, i.e., N(E = 0,T) ~ N(E,T = 0) with 
\E — Ep\ = kBT. This has been supported by Coulomb 
glass simulations^. We tested this relation by plotting 
N(E = 0,T) versus T, and N(E,T = T D ) versus E on 
the same graph, where T Q is the lowest temperature at 
which we were able to equilibrate the system. We show 
our results in Figure IT§1 for a — and 1. The hypothesis 
seems to work for a limited range of energies between 
ksT and the width of the Coulomb gap. It also appears 
to be more applicable for high disorder (er = 1) than for 
the case of no disorder (a = 0). 



E. Staggered Occupation 

We have studied the staggered occupation at various 
values of the disorder. At high temperatures the distri- 
bution has a peak centered at M s = for all values of the 
disorder. At low temperatures the distribution broadens 
and has two peaks symmetrically placed about zero for 
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FIG. 19: (a) N(E = 0,T) vs. T, and N(E, T — 0.120) vs. E 
for a = 0. The data is averaged over 45 runs, (b) N(E — 0, T) 
vs. T, and N(E,T = 0.0275) vs. E for a = 1. The data is 
averaged over 108 runs. The solid lines are guides to the eye. 



the ordered case and for small and moderate values of 
the disorder. For the strongly disordered case a > 0.5, 
the distribution has a peak centered at M s = for all 
values of the temperature where the system was able to 
attain equilibrium in our simulations. This is what one 
would expect for a random system. These features are 
illustrated in Figure I2D1 which shows the staggered occu- 
pation for various values of a in the vicinity of Tq- As 
a function of system size, the high temperature peak in 
P(M S ) becomes sharper as L increases for all values of 
a. An example is shown in Figure CHI 



V. SUMMARY 

We have performed a Monte Carlo study of a classical 
three dimensional Coulomb system of electrons in which 
we systematically increase the positional disorder by in- 
troducing deviations from positions in a cubic lattice. 
We start from a completely ordered system and gradu- 
ally transition to a Coulomb glass. The phase transition 
as a function of temperature is second order for all val- 
ues of disorder. We use finite size scaling to determine 
the transition temperature Tq and the critical exponent 
v. We find that Tc decreases and that v increases with 
increasing disorder. Both quantities saturate in the limit 
of large disorder. The specific heat peak value decreases 
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and the peak broadens to a broad bump with increas- 
ing disorder. A gap develops in the single particle den- 
sity of states for all values of a. At low temperatures 
N(E = 0) - T x where A > 3.8 for all values of a. At low 
temperatures and low energies near Ep, the density of 
states can be fit to a power law form N(E) <~ \E — Ep\ s 
where d — 1 < S < d for all values of a. S increases 
with increasing a, starting at S = 2.1 for a = and sat- 
urating at 5 — 2.65 for a = 1. The distribution of the 
staggered occupation has a single central peak at high 
temperature for all values of the disorder. In the ordered 
cases (tr < 0.4) P(M S ) develops two peaks symmetrically 
placed on either side of M s = in the vicinity of the 
phase transition. 



FIG. 20: Staggered occupation distribution for L = 8 for 
various values of a in the vicinity of To- o — (35 runs, 
T = 0.128), a = 0.1 (15 runs, T = 0.123), a = 0.2 (10 runs, 
T = 0.110), cr = 0.3 (10 runs, T = 0.085), a = 0.4 (40 runs, 
T = 0.045), a = 0.5 (10 runs, T = 0.03), and cr = 1 (10 runs, 
T = 0.03). The number of runs averaged over is indicated in 
parentheses. The solid lines are guides to the eye. 
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FIG. 21: Staggered occupation distribution for a = 1 at T = 1 
for L = 4, 6, and 8. The peak height increases with increasing 
L. The data shown is the result of averaging over 10 runs. 
The solid lines are guides to the eye. 
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